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Abstract 

A semi-numerical solution to Dokshitzer- Gribov-Lipatov- 
Altarelli-Parisi (DGLAP) evolution equations at leading 
order (LO), next-to- leading order (NLO) and next-to- next- 
to-leading order (NNLO) in the small- a; limit is presented. 
Here we have used Taylor series expansion method to solve 
the evolution equations and, t- and x-evolutions of the 
singlet structure functions have been obtained with such 
solution. We have also calculated t- and x-evolutions of 
deuteron structure functions F^, and the results are com- 
pared with the E665 data and NMC data. The results 
are also compared to those obtained by the fit to F$ pro- 
duced by the NNPDF collaboration based on the NMC 
and BCDMS data. 
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1 Introduction 

Structure functions in lepton-nucleon deep-inelastic scattering (DIS) are 
the established observables probing Quantum Chromodynamics (QCD), the 
theory of the strong interaction, and in particular the structure of the nu- 
cleon. The structure functions provide unique information about the deep 
structure of the hadrons and most importantly, they form the backbone of 
our knowledge of the parton densities, which are indispensable for analyses 
of hard scattering processes at proton- (anti-)proton colliders like the TEVA- 
TRON and the LHC. Structure functions are also among the quantities best 
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suited for precisely measuring the strong coupling constant ocg. The SLAC- 
MIT collaboration [1-3] observed the scaling behaviour [4] of the proton struc- 
ture function in DIS. This observation established the quark-parton model 
(QPM) as a valid framework for the interpretation of DIS data and these 
DIS processes can be expressed in terms of universal parton densities. In 
QCD, structure functions are defined as convolution of the universal parton 
momentum distributions inside the proton and coefficient functions, which 
contain information about the boson-parton interaction. At large momen- 
tum transfers Q 2 and not too small-x, where x is the fraction of proton 
momentum carried by the parton, QCD allows the perturbative calculation 
of the coefficient functions and predicts a logarithmic dependence (evolu- 
tion) of the proton structure functions with Q 2 to higher orders in as- Thus 
measurements of structure functions allow perturbative QCD to be precisely 
tested. The standard and the basic tools for theoretical investigation of DIS 
structure functions are the DGLAP evolution equations [5]. 

The solutions of the unpolarized DGLAP equation for the QCD evolu- 
tion of parton distribution functions have been discussed considerably over 
the past years [6-15]. There exist two main classes of approaches: those 
that solve the equation directly in x-space, and those that solve it for Mellin 
transforms of the parton densities and subsequently invert the transforms 
back to x-space. Some available programs that deal with DGLAP evolution 
are CANDIA [11] based on the logarithmic expansions in x-space, QCD PE- 
GASUS [12], which is based on the use of Mellin moments, HOPPET [14] 
and QCDNUM [15]. QCD PEGASUS which is a parton distribution func- 
tions [PDFs] evolution program based on Mellin-space inversion, has been 
used by the QCD Working Group to set some benchmark results [16-17]. 
HOPPET is an x-space evolution code that is novel both in terms of the 
accuracy and in terms of speed that it provides compared to other x-space 
codes. Among various methods used for the solution of this equation, most of 
the methods are numerical. Mellin moment space with subsequent inversion, 
Brute-Force method, Laguerre method etc. [18-20] are different methods 
used to solve DGLAP evolution equations. The shortcomings common to 
all are the computer time required and decreasing accuracy for x — > 0. For 
example, in the Brute-Force method the accuracy is better than 1% in the re- 
gion 10 -5 < x < 0.8. However, it takes significant amount of computing time 
[21]. The Laguerre method has an advantage of computing time. However, 
the accuracy is slightly worse at small x [22]. More precise approach is the 
matrix approach to the solution of the DGLAP evolution equations, yet it is 
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also a numerical solution. The matrix method is a matrix-based approach to 
numerical integration of the DGLAP evolution equations [23]. The method 
arises naturally on discretisation of the Bjorken x variable which is a neces- 
sary procedure for numerical integration. In Ref. [24], a quadratic contour 
has been derived for calculating the evolution of parton distribution functions 
within the Mellin transform method, and demonstrated its superiority over 
other techniques in the literature. With this method the parton distribution 
functions can be evolved according to any of the commonly used truncations 
of the evolution. Moreover, this method may also be used to extract par- 
ton distribution functions from collider data. Again, in Ref. [25], Stefan 
Weinzierl reports on a numerical program for the evolution of parton distri- 
butions which uses the Mellin-transform method with an optimized contour. 
An original implementation of a very elegant numerical treatment used to 
solve the DGLAP evolution equations at NLO is presented in Ref. [26]. This 
method consists in expanding parton distributions and splitting functions on 
Laguerre polynomials, which reduces DGLAP equations to a set of ordinary 
differential equations defined by recurrence. Recently, the NNPDF collabo- 
ration [27] presented and utilized an improved fast algorithm for the solution 
of evolution equations and the computation of general hadronic processes. 
As an alternative to the numerical solution, one can study the behavior of 
quarks and gluons via analytic solutions of the evolution equations. Although 
exact analytic solutions of the DGLAP equations cannot be obtained in the 
entire range of x and Q 2 , such solutions are possible under certain conditions 
[28-29]. Some approximated analytical solutions of DGLAP evolution equa- 
tions suitable at small-x, have been reported in recent years [6-9] [30-32] with 
considerable phenomenological success. For example, in Ref. [6] the spin- 
dependent singlet and non-singlet structure functions have been obtained by 
solving DGLAP evolution equations at LO and NLO in the small-x limit by 
using Taylor series expansion and then the method of characteristics. They 
have calculated the t- and ^-evolutions of deuteron structure function and 
the results show very good agreement with the SLAC E-143 Collaboration 
data. Similarly, in Ref. [7], the spin- independent singlet and non-singlet 
structure functions have been obtained by solving the DGLAP evolution 
equations by using Taylor expansion and then the method of characteristics 
at LO and NLO. Here also the computed values of the t- and re-evolutions of 
deuteron structure function are in good agreement with the NMC data. An 
approximate solution of the DGLAP equation by using Taylor expansion is 
presented and the x-distribution of the deuteron structure function is calcu- 
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lated at small- x in Ref. [10]. Here The results are compared with the EMC 
NA 28 experiment data and it shows good phenomenological success. The 
one-loop splitting functions corresponding to LO DGLAP equation is given 
in Ref. [33]. Similarly the two- loop splitting functions governing the evo- 
lution have been known for a long time [34]. The calculation of the NNLO 
QCD approximation for the structure functions F 2 of DIS is important for 
the understanding of perturbative QCD and for an accurate comparison of 
perturbative QCD with experiment. To obtain the NNLO approximation for 
these structure functions one needs the corresponding 3-loop splitting func- 
tions. Recently the three-loop splitting functions governing the evolutions 
of both the non-singlet and singlet structure functions have been reported 
which are in with good agreement with the existing data [35-40]. 

In a recent paper [8] the non-singlet structure functions have been ob- 
tained by solving DGLAP evolution equations at NNLO in the small- x limit. 
Here a Taylor series expansion has been used to solve the evolution equations 
and to obtain the semi-numerical solution. The results are compared with 
the Fermilab experiment E665 data and New Muon Collaboration (NMC) 
data. In the present paper we extend the work to the singlet case. Here 
we intend to solve the DGLAP equation analytically by Taylor expansion 
method to obtain a semi-numerical solution for singlet structure functions at 
small- £ at NNLO. In this work we have also calculated t- and ^-evolutions 
of deuteron structure functions and the results are compared with E665 data 
[41] and NMC data [42]. We have also compared our results with those ob- 
tained by the fit to F 2 produced by the NNPDF collaboration [27, 43] in 
2002 based on the NMC and BCDMS data. The NNPDF parametrizations 
have presented a determination of the probability density in the space of 
structure functions for the F 2 structure function for proton, deuteron and 
non-singlet, as determined from experimental data of the NMC and BCDMS 
collaborations. Their results take the form of a set of 1000 neural nets, for 
each of the three structure functions, each of which gives a determination 
of F 2 for given x and Q 2 . The central value and the errors of the structure 
functions determined in the NNPDF fit can be computed out of the ensemble 
of 1000 nets according to standard Monte Carlo techniques. 

The spin-dependent structure functions of the nucleon, g\ and g 2 , can be 
measured in polarized DIS [44] by using polarized lepton beams scattered by 
polarized targets. The first experiment in polarized electron-polarized proton 
scattering was performed in the 1970s which helped to establish the parton 
structure of the proton. Polarized DIS lepton nucleon scattering experiments 
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have been performed mainly at CERN, SLAC, DESY, and JLAB, and these 
processes have played a key role in understanding QCD and the spin structure 
of the nucleon [45-51]. The solutions of spin-dependent DGLAP evolution 
equations give polarized quark and gluon parton densities which ultimately 
give polarized proton and neutron structure functions. The Taylor series 
expansion method that gives acceptable solutions for spin-independent evo- 
lution equations [8,10], which we use in our present work, can also be used 
for the spin-dependent case [52]. The LO evolution equation of polarized 
structure functions implies scaling violations. The NLO and higher order 
evolution equations of polarized structure functions will help us to under- 
stand the higher twist contributions to the scaling violations and, ultimately, 
the spin content of the nucleon. 



2 Solution of singlet structure function at 
NNLO 

The singlet quark density of a hadron is given by [36] 

N f 

q s (x, /i/ 2 , n r 2 ) = ^2[qi(x, /i/ 2 , n 2 ) + ft(ar, /if 2 , /i r 2 )], (1) 
i=i 

where qi(x, fif 2 , /i r 2 ) and qi(x, fif 2 , fi 2 ) represent the number distribution of 
quarks and antiquarks, respectively, in the fractional hadron momentum x. 
The corresponding gluon distribution is denoted by g(x, fif 2 , fi 2 )- The sub- 
scripts i indicate the flavor of the quarks (and antiquarks), and Nf is the 
number of effectively massless flavours, and finally \x T and fif represent the 
renormalization and factorization scales respectively. Since fi r and fif are 
arbitrary, their values may be chosen independently. 

The DGLAP evolution equation in the singlet sector in the standard 
form is given by 



9 ( 9s \ = ( P m P q9 \(cx( qs 
dlnQi \ g J \P gq P gg J^y\ g 

where (^) represents the standard Mellin Convolution in the momentum vari- 
able and the notation is given by 

a(x)(g)b(x)= f l ^a(u)b(-) (3) 
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Thus, using equation (3), equation (2) can be written as 

d ( qs\ f'duf P qq P qg \fqs 
dlnQ* \ g J J x uj \P gq P ag )\ g 

(4) 

Hi> i * ggV / 

3qs f^(P, tqs + P v9 ) (5) 

Jx W 



1 du ( P qq q s + P qg g 

UJ V P 9q1S + Pgg9 



which implies, 



<91nQ 
dg f 1 du 



[x 



dlnQ 2 J/i( P ^ + P ^ («) 

where, qs and g are singlet and gluon parton densities respectively and P qq , 
P qg , Pgq-, P gg are splitting functions. The quark-quark splitting function P qq 
can be expressed as [53] 

+0(P<*\x)) (7) 

where, P qq \x), P q l\x) and P qq \x) are LO, NLO and NNLO quark-quark 
splitting functions respectively. Other splitting functions can be expressed 
in the same way. 

The splitting functions can be obtained from the N-space results of the 
Mellin space by an inverse Mellin transformation. The quark-quark splitting 
function P qq in equation (2) can be expressed as P qq = P^ s + Nf(P qq + P^) = 
P£ s +Pps- Here P^ s is the non-singlet splitting function calculated up to the 
third order [38]. P qq and P^ q are the flavor- independent contributions to the 
quark-quark and quark-antiquark splitting functions respectively. At very 
small- £ the pure singlet term Pp $ dominates over P^ s [40] . The gluon-quark 
(Pgq) and quark-gluon (P qg ) entries in equation (2) are given by P qg = NfP qi9 
and Pgq = P gqi , where P qig and P gqi are the flavor- independent splitting 
functions. 

Substituting the LO splitting functions in equation (5) and simplifying, 
the DGLAP evolution equations for singlet structure function at LO can be 
written as 

^^-^g{3 + 41n(l-a;)}Ff(x,t) + /f(x,t)] =0, (8) 
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where, 



N r 



Ff(x, t) = q s = + 
4 Z" 1 



(9) 



i=i 



lf(M) = | / -^-[(1 + u, 2 )if (-,t) - 2if(x,t)] 

O J x 1 — U! 

+N f J {u 2 + {l-u) 2 }G(^-^du. 



(10) 



Q 2 

A 2 



Here i = In— ,A is the QCD cut off parameter. It is the scale at which 



partons turn themselves into hadrons. Also G(^,t \ = g is the gluon parton 
densities. 

Similarly, substituting the splitting functions upto NLO and upto NNLO 
in equation (5) and simplifying, we get the DGLAP evolution equations for 
singlet structure function at NLO and NNLO respectively as, 



dF£(x,t) a s (t) 



Of 



2tt 



-{3 + 4 ln(l - x)}F 2 s (x, t) + if (x, t) 



2vr 



l 2 5 (x,t) = 0, 



(11) 



dFiM ^ (t) [|{3 + 41n(l-x)}Pf(x,t) + if(x,t) 



Of 



2tt 



-(^M-(^)Wo 



;i2) 



where, 



I*(x,t) = (x-l)FZ(x,t) f(u)du>+ f(u)F£(—,t)(L} 

J0 Jx u 

'^(^(^j^ + ^^MG^)^ (13) 



lf(x,t) 



1 do; 
a; 



Pqq( x )F 2 



t)+P qg (x)Gl-,t 



X 



(14) 



The explicit forms of the third order splitting functions are given in 
Appendix A. 
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The strong coupling constant, otsiQ 2 ) is related with the /3-function as 

[54] 

^ = "d%Q^ = - 16^ - 64^ + ° M (15) 

where, 

fo = T Ne ~t T, = u ~l Nf ' (16) 

34 1 D 38 

ft = -jA? - jN c N, - 2C F N, = 102 - —AT, (17) 

ft = 2 -^N? + 2C%T, - f Cr N e T, + fc F T} + ^§N^ 

= 2857 _ 6673 325 

6 18 ' 54 y v ; 

are the one-loop, two-loop and three-loop corrections to the QCD /3-function 

N 2 -l 

and Nf being the number of quark flavor. We have used N c = 3, Cp = — — 

= - and Tt = -Nf. 

3 1 2 1 

Now let us introduce the variable u = 1 — oj as discussed earlier [6-8]. 
Therefore, we can write ^ = (jz^y = x + jz^- Thus we can expand F^ix^t) 
by Taylor expansion method as 



F 2 s (- 7 t) = F 2 s (x + -^-,t) 
00 1 — u 



xu \dFn(x,t) 1 / xu \ 2 d 2 Fn(x,t) , , 

Fi(x,t) + [- V — + nh — isM J - + --- ( 19 ) 



it/ <9x 2 Vl — u/ 9 2 x 

As x is small in our region of discussion, the terms containing x 2 and 
higher powers of x can be neglected as those terms are still smaller and 
therefore, we can rewrite 

Similarly, 

G( £, i) = G (*, t ) + -^M. (2i) 

1 — "U OX 
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Using equations (20) and (21) in equations (10), (13) and (14), and 
performing w-integrations, equations (8), (11) and (12) become 



dF 2 s (x,t) as (t) 



dt 



2tt 



A 1 {x)F 2 s {x,t) 



x) dF i {x ' t] +A 3 (x)G(x,t) 



dx 



+AJx 



dG(x,t)- 
dx 



(22) 



dF 2 s (x,t) a s (t) 



dt 



2n 



+AJx) 



A^Fffat) + A 2 (x) dF ^ X,t) + A 3 (x)G(x,t) 



dG(x,t)i t , as (t) 



-B»(x 



dx 



dFj(x,t) 
dx 



2tt 



BAx 



dx 



Bi{x)Fi{x,t) 



dG(x,t) 
dx 



(23) 



dF 2 s (x,t) a s (t) 



dt 

+AJx] 



2tx 
dG(x,t) 
dx 



A 1 (x)F 2 s (x,t) + A 2 {x) dF ^ X ; t] + A 3 (x)G(x,t] 



a s (t) 



2ix 



dx 

B 1 (x)F 2 s (x,t)+B 2 (x 



+B 3 (x)G(x,t) + B A {x 
dF 2 s (x,t) 



dG{x,t) 
dx 



+ 



as(t) \ 3 
2ix J 



+C 2 (x)- 



dx 



+ C 3 (x)G(x,t) + C i (x) 



dF£(x,t) 
dx 

C 1 (x)Fi(x,t) 
dG(x,t) 



dx 



(24) 



where, A^x), Bi(x) and Ci(x) (where 2=1,2,3,4) are functions of x (see 
Appendix B). 

The Q 2 -evolution of the proton structure function F 2 (x, Q 2 ) is related to 
the gluon parton densities in the proton, G(x, Q 2 ), and to the strong interac- 
tion coupling constant, as- The gluon parton densities cannot be measured 
directly through experiments. It is, therefore, important to measure the 
gluon parton densities G(x, Q 2 ) indirectly using the proton structure func- 
tions F 2 (x,Q 2 ). Hence the direct relations between F 2 (x,Q 2 ) and the gluon 
parton densities G(x, Q 2 ) are extremely important; using those relations the 
experimental values of G(x, Q 2 ) can be extracted using the data on F 2 (x, Q 2 ). 
Therefore, in the analytical solutions of DGLAP evolution equations for sin- 
glet structure functions or gluon parton densities, a relation between singlet 
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structure function and gluon parton densities has to be assumed [6,7]. The 
commonly used relation is G(x,t) = K{x)F 2 {x, t), where K(x) is an ad hoc 
function of x. Since these evolution equations of gluon parton densities and 
singlet structure functions are in the same forms of derivative with respect 
to t, so we can consider this function. And also the input singlet and gluon 
parameterizations, taken from global analysis of parton distribution func- 
tions, in particular from the MSTW08 parton set, to incorporate different 
high precision data, are also functions of x at fixed Q 2 [55]. So the rela- 
tion between singlet structure function and gluon parton densities will come 
out in terms of x at fixed-Q 2 . However the actual functional form of K(x) 
can be determined by simultaneous solutions of coupled equations of singlet 
structure functions and gluon parton densities. 

Hence equations (22), (23) and (24) take the form 

_ t dFipt) + Li(M) AF£M + Ml{Xjt)F s M = q 5 (25) 

dF s (r t) 

L 2 (x,t) 1 + M 2 (x,t)F 2 s (x,t) = 0, (26) 

L 3 ( x , t) dF "^ t] + M 3 (x, t)F 2 s (x, t) = 0, (27) 
where 

L x {x,t) = ^A f (t)[A 2 (x) + K(x)A 4 (x)}, (28) 

M 1 (x,t) = ^A f (t)[A 1 (x) + K(x)A 3 (x) + ^i 4 (.)], (29) 

L 2 (x,t) = ^A f (t)[(A 2 (x) + K(x)A A (x)) + T (B 2 (x) + K(x)B A (x))}, (30) 

M 2 (x,t) = ^A f (t)[{A 1 (x)+K(x)A 3 (x) + ^A 4 (z)) + T Q (B 1 (x) 

+K(x)B 3 (x) + ^^-B4x))}, (31) 

L 3 (x, t) = ^A f (t)[(A 2 (x) + K{x)A i {x)) + T (B 2 (x) + K(x)B 4 (x)) 

+T 1 {C 2 (x) + K{x)C 4 (x))}, (32) 

M 3 (x,t) = ^A f (t)[{A 1 (x)+K(x)A 3 (x) + ^^A 4 {x)) + T Q (B 1 (x) 



<9if (x, 


t) 


dt 




dFi(x, 


t) 


dt 




dFi(x, 


t) 


dt 
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+K(x)B 3 (x) + f^Mfl 4 (x)) + T^dix) + K(x)C 3 (x) 



and 



dx 



a s (t) _ 3A f (t) 



(33) 



(34) 



2ir It 

Also one can consider two numerical parameters To and T\, such that 
T 2 (t) = T T(t) and T 3 (t) = T{T{t), where T(t) = As discussed in [8], 

here also we have considered To = 0.05 and T\ = 0.006 within the range 
< Q 2 < 30 GeV 2 . These values are chosen such that differences between 
T 2 (t) and T T(t), and also T 3 (t) and T\T(t) are negligible. Thus the consid- 
eration of the parameters T an d Ti does not give any abrupt change in our 
result. 

The general solution of the equations (25), (26) and (27) is F(U, V) = 0, 
where F(U, V) is an arbitrary function. Here, U (x, t, i 7 ^) = K\ and V(x, t, T 2 5 ) 
= K 2 are two independent solutions of the equation 



dx 



dt 



dFi ( 



x,t) 



Li(x,t) -t -Mi(x,t)Fi(x,t) 



(35) 



To get the solution of equations (25), (26) and (27) we will follow the same 
procedure as discussed in Ref. [8]. Now, let us introduce functions Li(x) and 
Wi(x) such that, U{x,t) = &T(t)Ll(x) and Mi(x,t) = ^T(t)M~i(x) where 
i=l,2,3 for LO, NLO and NNLO respectively. 

Solving this equation for LO, we obtain 



U(x, t, F 2 S ) = t ■ exp 



Li(x) 



-dx 



(36) 



and 



where 



V(x,t,F£) =F£(x,t)-exp 



ill ■ 



L 1 {x) = — (A 2 + KA 4 ), 
Po 

W 1 (x,t) = j o (A 1 + KA 3 + ^A 4 ). 



(37) 

(38) 
(39) 
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Thus it has no unique solution. The simplest possibility to get a solu- 
tion of the equation (25) is a ■ U + /3 ■ V = 0, where a and (3 are arbitrary 
constants. Putting the values of U and V in this equation we get 



at ■ exp 
from which we get 



— dx 
LAx) 



+ f$F§(x,t) ■ exp 



LAx) 



0. 



(40) 



F 2 s (x,t) = -Tt-expT / ( 



MAx) 



LAx) LAx) 



)dx 



where 7 = % is another constant. Now defining 



(41) 



F 2 s (x,t ) 



-7*0 • exp 



K LAx) LAx)' ■ 



(42) 



Q 2 

at t = t , where t = In — , for any lower value Q 2 = we get from equa- 
tion (41) 

(43) 



t 



F 2 s (x,t) = F 2 s (x,t )(-) 



Again defining at x = xq, 



F 2 s (x ,t) = --ft -exp / ( 



MAx) 



l Li(x) Li(x) 



)dx 



X=XQ 



(44) 



we obtain 



F^x.t) = F 2 s (x ,t) - exp 



xo Li(x) L x (x) 



(45) 



Equations (43) and (45) give the t- and x- evolutions of singlet structure 
functions at LO, respectively. 

The deuteron structure function can be written in terms of the singlet 
quark distribution function as [3], 



F 2 d (x,t) = °-F 2 s (x,t) 



(46) 
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Substituting equations (43) and (45) in equation (46), the t- and x- 
evolutions of deuteron structure function at LO can be obtained as 



and 



F«(x,t) = F 2 d (x,t )(l) 
to 



F d (x,t) = F d (x A)-exp / ( 



MAx) 



)dx 



(47) 



(48) 



, Xo 'Li(x) Li(x) 

Proceeding in the same way we obtain the t and x-evolution of deuteron 
structure function F d (x, t) at NLO as 



F d (x,t) = F d (x,t ) 



fl+b/t 



1+V*0 




■ exp 



fb b 

\t~T 



and 



where 



F 2 d (x,t) = F 2 d (x A)-e W I (l-MM^x 

L 2 [X) L 2 {X) 



X,) 



(49) 



(50) 



li-^-^Mjijk-W^ (51) 



f (l+6/t) 



exp [ - J ■ exp 



L,(x) 



L 2 (x) \x=xq) 



L 2 (x, t) = —[(A 2 + KA A ) + T (B 2 + KB A )}, 

Po 



(52) 
(53) 



O r-) 7, ^ O TV*" 

W 2 (x,t) = j[(A x + KA 3 + —A,) + To^x + KB 3 + ^5 4 )], (54) 

Similarly, the t- and x-evolutions of deuteron structure function at NNLO 
can be obtained as 



t i+{ b ~b 2 )/t 

F d (x,t) = F d (x,to){ {b _ b2)/to 

^2 1„2 



• exp 





b-c-tf\x\ t 



b-c-b 2 \n 2 t 



to 
13 



(55) 



and 



where 



F*(x,t) = F 2 d (x ,t)-exp 



L 3 (x) L 3 (x) 



(56) 



.(l+(6-6»)/*o) 
7 r o 



cxp 



b - c - b 2 In 2 tr- 



io 



)- exp [/ ( i: 



L 3 (x) 



■&^)dx]} 



(57) 



F 2 d (x ,t) 



5|_ 7t (i + (^)A). exp ^---& 2 ln * 



■ exp 



L 3 (x, t) = —[{A 2 + KM) + T (5 2 + KB 4 ) + Ti(C 2 + ATC 4 )], 



A) 

M^x, t) = + -ft' A3 + ^-A 4 ) + To (fix 

8 K 

+T 1 {C 1 + KC 3 + —C 4 )}, 



(58) 
(59) 



KB, 



dK 

dx 



B A ) 



(60) 



with 6= 4 and c= §• 

Here the input for F£(x, to) is taken from experimental data correspond- 
ing to the lowest-Q 2 point. Similarly, the input for F$(xQ,t) is taken from 
the experimental data corresponding to the highest- x point. 

HERA shows that the structure function has a steep behavior in the 
small-x region 10~ 5 < x < 10~ 2 . This steep behavior is well described in 
the framework of the DGLAP evolution equations. So we can expect that 
the method we used in this work to solve the DGLAP evolution equations 
is valid in the small-x region; roughly in the region 10~ 3 < x < 10 _1 . This 
method may also be applied in the region where x < 10 -3 , but we have not 
checked the validity of the method in this region. However, our method 
is not valid at very small-x where recombination processes between gluons 
have to be taken into account, since DGLAP equations fail to describe the 
recombination processes. 
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3 Result and discussion 



In this paper we have calculated the t- and x-evolutions of singlet struc- 
ture functions. The deuteron structure function is directly related to the 
singlet structure function as given by the equation (46). To determine the 
proton structure function we need to know both the singlet and non- singlet 
structure functions. Therefore in this particular work we have compared the 
results only to the deuteron structure function measured in fixed-target ex- 
periments. We have compared our results of t- and re-evolutions of deuteron 
structure function with NMC and E665 experimental data in figure l(a,b) 
and figure 2(a,b) respectively. In figure 3(a,b) we have compared our results 
with thosed obtained by NNPDF parameterizations. We consider the range 
0.0052 < x < 0.18 and 1.094 < Q 2 < 26 GeV 2 for E665 data, 0.0045 < x 
< 0.18 and 0.75 < Q 2 < 27 GeV 2 for NMC data. We have used the range 
0.0045 < x < 0.09 and 1.25 < Q 2 < 26 GeV 2 to compare our results with 
those obtained by the fit to F£ produced by the NNPDF collaboration in 
2002 based on the NMC and BCDMS data. In the t-evolutions of deuteron 
structure function our computed values of F^{x, t) from equations (47), (49) 
and (55) for LO, NLO and NNLO respectively are plotted against Q 2 for 
different values of x. On the other hand, for x-evolutions our computed val- 
ues of F$(x,t) from equations (48), (50) and (56) for LO, NLO and NNLO 
respectively are plotted against x for different values of Q 2 . For convenience, 
value of each data point is increased by adding 0.5* or OAi or 0.3* or 0.2i or 
O.li where % = 0, 1, 2, 3, • • • are the numberings of curves counting from the 
bottom of the lowermost curve as the 0-th order. In all graphs, the lowest-Q 2 
and highest-x points are taken as inputs for F 2 d (x, t Q ) and F 2 d (x ,t) respec- 
tively. Here vertical error bars are total statistical and systematic uncertain- 
ties. It is observed that the best-fit curves are obtained for 2 < K < 8. Here 
we have also performed a chi-square test of our computed results and we got 
chi-square per degree of freedom of about 0.9. 

We observe that the t- and x-evolutions of deuteron structure function 
are in good consistency with the experimental data and parametrizations. As 
observed from this work it is clear that the region of validity of our method 
is approximately in the range 10 -3 < x < 10 _1 and 0.5 < Q 2 < 30 GeV 2 . 
But this methed is also valid for other ranges of x and Q 2 . Though various 
methods like Laguerre polynomials, Brute-Force method, Mellin transforma- 
tion etc. are available in order to obtain a numerical solution of DGLAP 
evolution equations, our method to solve these equations analytically is also 
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a workable alternative. Here we consider a parameter K(x) in assuming a 
relation between singlet structiure function and gluon parton densities. We 
have also used two other parameters like To and T\. However the number 
of parameters used here is smaller compared to other methods where several 
parameters have been used mainly in input functions [55, 56]. Moreover with 
this method we can calculate the x-evolutions of deuteron structure function 
in addition to the t-evolutions. 
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Appendix:A 

In Ref. [40] the three loop splitting functions have been presented in 
both Mellin-N and Bjorken-x space. The exact expressions given in equations 

( 2 1 ( 2 1 

(4.12) and (4.13) in reference [40] for the functions P P g(x) and P\ g (x) are 
not simple to use. Therefore they have also presented compact approximate 
representations. As it is shown in reference [40], at small-x (x — > 0) the 

(2} (2\ 

coefficients of the 1/x terms in the definition of Pps(%) and Pq 9 (x) can be 
written in terms of the color factors and the Riemann ^-functions. After 
inserting the values of the color factors and the approximated numerical 

(2) / (2) 

values of ^-functions, the three loop splitting functions P p ^(x) and Pq 9 (x) 
are approximated in the limit (x — > 0) as given in Ref. [40]. The third-order 
pure-singlet contribution to the quark-quark splitting function is 



f£2(*) 



N f (-5.921% - 9.751L 2 - 72.11Li + 177.4 + 392.9x - 101.4x 2 



-57.04L L! - 661.61L + 131.4L 2 - — L% + - 506.0x _1 

9 27 

■^ir-%) + N 2 f (1.778Ll + 5.9441^ + 100.1 - 125.2a; + 49.26x 2 - 12.59x 3 
27 

32 256 i 

-1.889L Li + 61.75L + 17.89L 2 + — L 3 + —x^ 1 ) (1 - x), (61) 

27 81 J 
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The three-loop quark-gluon splitting function is 



100 



70 



P$(x) 2* N f (-^-L\ - -^-L\ - 120.5L^ + 104.42^! + 2522 - 3316x + 2126a;' 



44 

-L L!(1823 - 25.22L ) - 252.5xL 3 + 424.9L + 881.5L 2 , - — L 3 



536 

'^7 



L{- 1268.3a;- 1 - 5^"%) +JV?(^L? 



200 

"27" 



L\ - 5.496L! 



-252.0 + 158.0x + 145.4x 2 - 98.07xL 2 Q + 11.70xLg - L £i (53.09 



+80.616L ) - 254.0Lo - 90.80L; 



376 



-ljr 



9 



243 



(62) 



The non-singlet splitting function calculated upto third order is given by 



N f 



{L^-lGSSx- 1 - 7.208x) + 151.49 + 44.51x - 43.12x' 



+4.82x 3 }(l -x) + LoLi(-173.1 + 46.18L ) + 178.04L 

40 

+6.892L 2 + -(L A - 2Ll) 
Again, the functions F g (co) and F g (co) are defined as 



(63) 



F ? » = 2C F T R N f F qq (co), 

- C F T R N f FUu>) + C G T R N f F 2 (uj) 



where, 



(64) 
(65) 



20 



56 



2 + 6uj- -co 2 

9co 9 

-(1 + co)ln 2 (co) 



/ 8 \ 

(l + 5co + -co 2 ) ln(co) 



FUuj) = 4 - 9co - (1 - 4w) ln(w) - (1 - 2w) ln 2 (u;) + 41n(l - to) 



21n 2 (- 

to 

182 14 
"9" 4 



CO , 



41n( 



40 

— co H 

9 9a; 



,136 



-7T 

3 

38 



10 



-a; 



)ln(w) -41n(l - w) 



44 



-(2 + 8w) ln 2 (w) + - ln 2 (w) + — ln(w) - 2 ln 2 (l - w) 

o 



(66) 



(67) 
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+41n(l-w) 



7T 2 218 



" z z 



(68) 



Also, the function f(co) is given as 



2N r 



where, the two-loop quark-gluon splitting function Pqg is defined in Ref. 
[34]. Here, T R = |, C G = N c = 3 and C F = = |. 



/(w) = Cf.[P F (a;) - + -C F C A [P G + P A (w)]+ 

C F T R N f P Nf (u;), 



where, 



PnAu) 



i±^(-lnc-|)-2(l-c) 

1 — Cu o 



(69) 
(70) 



2(1 + u 2 ) 
' (1- W ) 



ln(w) ln(l - u) 



1-LO 



2lo ) In uj 



1 40 
--(l + w)lnu;+y(l-u;) 



(71) 



J-* / \ (1 +w 2 ) /, 2 , , 11 w x 67 7T 2 \ 1 ^ 

P G ( W ) = ( ln a ( W ) + T In( W ) + - - T ) - -(1 + W ) Inc. 



+ T (l-«) 



(72) 



2(1 +u 2 ) dk. n-k\ 

P A (u) = -7r-, — r- I — 111 (— — ) + 2(l + w)ln(w) 



(1 + w) 



k 

-4(1- w) 



k 



(73) 

with Lq = ln(x), L\ = ln(l — x). Here results are from direct x-space 
evolution and are calculated using Fortran package [35-40]. 

Appendix:B 

The explicit forms of the functions Ai(x), Bi(x) and Ci(x) (where z=l,2,3,4) 

are 



Ai(x) = 2x + x 2 + 41n(l - x), 



(74) 
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A 2 (x) = x — x 3 — 2xha.(x), 

A 3 (x) = 2N f (^-x + x 2 -^x 3 ), 

5 2 
AAx) =2N f (--x + 3x 2 - 2x 3 + -x 4 - x \n(x)) . 

v 3 3 ' 

4 



i(x) = x f(u)du- f(u)du + -N f F qq (u)du 

JO JO ^ Jx 

B 2 {x) =xf [f(u) + UfF^cu)] —du, 

Jx 6 U 

B s (x)= f F q s g (u)du, 

J X 
f 1 

B 4 ( x ) = x / -Z^-F q s g {u)(kj, 

Jx U 

C 1 {x) = N f ——R^), 

Jo L—OJ 

^) = N f / ^(w), 

Jo l 1_w J 

ujxduj 

o (T 3 ^ 



Z 1 x uxdu . , 
^) = / — — -^R 2 {u) 
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where, 



R^u) = {ln(w) ln(l - to) [-173.1 + 46.18 ln(l - u>)] + 178.04 ln(l - u) 
+6.892 ln 2 (l -u) + §[hi 4 (l - w) - 2 ln 3 (l - w)]} + w{ln(w) 
(-163.9(1 - a;)" 1 - 7.208(1 - to)) + 151.49 + 44.51(1 - u) 
-43.12(1 - to) 2 + 4.82(1 - to) 3 } + w 2 {-5.926 ln 3 (u;) 
-9.751 \n 2 (tu) - 72.11 ln(w) + 177.4 + 392.9(1 - to) 
-101.4(1 - uo) 2 - 57.041n(l - w) ln(w) - 661.6 ln(l - w) 
+ 131.4 ln 2 (l - oj) - *f ln 3 (l - u) + *f ln 4 (l - w) 
-506.0(1 - u/T 1 - ^p(l - w)" 1 ln(l - w)} + A^l.778 ln 2 (w) 
+5.944 ln(w) + 100.1 - 125.2(1 - u) + 49.26(1 - w) 2 
-12.59(1 - w) 3 - 1.889 ln(l - u>) ln(w) + 61.75 ln(l - w) 

+ 17.89 ln 2 (l + g ln 3 (l - w) + fp(l - w)- 1 } (86) 

i22( w ) = { ^ ln 4 (w) - f ln 3 (u;) - 120.5 ln 2 (w) + 104.42 ln(w) + 2522 
-3316(1 - w) + 2126(1 - w) 2 - 252.5(1 - w) ln 3 (l - to) 
+ ln(w) ln(l - uj) (1823 - 25.22 ln(l - uo)) + 424.9 ln(l - u) 
+881.5 ln 2 (l - u) - f ln 3 (l - w) + ^ ln 4 (l - w) - 1268.3 
(1 _ w )-i _ 8|6(i _ w )-i l n (i - w )} + iV f {§ ln 3 (u;) + 202 ln 2 (w) 
-5.496 In(w) - 252.0 + 158.0(1 -u) + 145.4(1 - w) 2 
-139.28(1 - cu) 3 - 98.07(1 - to) ln 2 (l - w) + 11.70(1 - u>) 
ln 3 (l - u) - 1b(w) ln(l - w)(53.09 + 80.616 ln(l - w)) 
-254.0 ln(l - w) - 90.80 ln 2 (l - cu) - ^ ln 3 (l - u) 

-f ln 4 (l-c) + ^(l-u;r 4 } (87) 
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Figure 1: NMC data. Dotted lines are LO results, dashed lines are NLO and 
solid lines are our NNLO results. For clarity, data are scaled up by +0.2i (in 
Fig. (a)) and +i (in Fig.(b)) (with i = 0,1,2,3) starting from the bottom of 
all graphs in each figure. 
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Figure 2: E665 data. Dotted lines are LO results, dashed lines are NLO and 
solid lines are our NNLO results. For clarity, data are scaled up by +0.4i (in 
Fig. (a)) and +0. li (in Fig.(b)) (with i = 0,1,2,3) starting from the bottom 
of all graphs in each figure. 
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Figure 3: NNPDF data. Dotted lines are LO results, dashed lines are NLO 
and solid lines are our NNLO results. For clarity, data are scaled up by +0.4i 
(in Fig.(a)) and +0.1i (in Fig.(b)) (with i = 0,1,2,3) starting from the bottom 
of all graphs in each figure. 
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